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Abstract 

We formulate a random matrix model which mimics the chiral phase transition in QCD 
with two light flavors. Two critical exponents are calculated. We obtain the mean field 
values /3 = I and 6 = 3. We also find that the chiral phase transition can be characterized 
by the dynamics of the smallest eigenvalue of the Dirac operator. This suggests an 
alternative order parameter which may be of relevance for lattice QCD simulations. 
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1 Introduction 



In recent years the QCD phase transition has been studied in a variety of ways using both 
numerical simulations and analytical methods. Such studies have led to the conviction 
that it may be a second-order phase transition [|l], ^. As has been stressed, in particular in 
P, ^ , this has important consequences because the transition can then be characterized by 
critical exponents corresponding to a specific universality class. In particular it was argued 
that the critical exponents are those of an 0(4) Heisenberg spin model. (On a lattice and 
with Kogut-Susskind fermions, the relevant model might rather be an 0(2) spin model.) 
However, as always, universality arguments must be used with care. According to a recent 
suggestion in [^j based on simulations of the three dimensional Gross-Neveu model, they 
may not be valid for phase transitions involving soft modes composed of fermions. The 
reason is that the lowest Matsubara frequency suppresses infrared divergencies which 
lead to the universal critical exponents. Indeed, Kocic and Kogut found that the critical 
exponents in their model are given by mean field theory. 

We wish to study the chiral phase transition from the perspective of the spectrum of 
the Euclidean Dirac operator. Although initial numerical lattice results regarding this 
issue have become available [P, ^, 0, a systematic study is beyond the reach of present 
day computers. Thus, we prefer to construct a simple model which contains the main 
ingredients of QCD related to chiral symmetry breaking. The zero temperature version of 
this model was considered previously [§, and it was shown that it is equivalent to the 
finite volume effective partition function. In particular, it was shown that the spectrum 
of the Dirac operator in this random matrix model obeys the so-called Leutwyler-Smilga 



sum rules 10 



In this work we consider the non-zero temperature results of two versions of the random 
matrix model, for SU{2) and for SU{Nc), Nc > 3 both with fundamental fermions, 
which are characterized by real and complex matrix elements, respectively. We shall 
give analytical results for the critical exponents j3 and 6 as well as numerical results (for 
the case of complex matrices only). We also study the temperature dependence of the 



spectrum. Finally, we calculate the dependence of the chiral condensate on the valence 
quark mass. These results will be shown to be in qualitative agreement with recent lattice 
calculations [|2|. 

2 Finite temperature chiral random matrix model 

The QCD partition function for vacuum angle 6 is defined as 

ZQCD = Y.^''''ZQCD{y) , (2.1) 

V 

where the partition function in a sector with topological charge v and A^/ fermionic flavors 
is given by 

Nf 

ZQCD{iy) = {Y[^} 11(^1 + ^})) A ■ (2.2) 

/=1 A„>0 

Here, (■ ■ ■)a denotes the averaging over the gauge field configurations, with topological 
charge u, weighted according the the QCD action. The eigenvalues, A^, of the Dirac op- 
erator fluctuate over the ensemble of gauge field configurations. In general, the complete 
eigenvalue density is a determined in a nontrivial way by the QCD dynamics. How- 
ever, fluctuations of the eigenvalues on the microscopic scale (on the scale of the average 
level spacing) are believed to show universal characteristics. It is our conjecture that the 
eigenvalues near zero virtuality obey such microscopic universality. This implies that the 
detailed dynamics of the QCD partition function are not important for such fluctuations, 
and can be described equally well with a random matrix ensemble which respects the 
global symmetries of the Dirac operator. In particular, the following properties are in- 
cluded IQ: (i) The Ua^^) symmetry leading to an eigenvalue spectrum ±A„. (ii) The zero 
mode structure of the Dirac operator. (In the sector of topological charge u the Dirac 
operator has exactly u zero eigenvalues all of the same chirality.) (iii) The flavor chiral 
symmetry and its spontaneous or explicit breaking, (iv) The reality type of the represen- 
tation of the gauge group. For SU{Nc), Nc > 3 in the fundamental representation, the 
gauge field is complex and so are the matrix elements of the Dirac operator. The gauge 
group SU (2) in the fundamental representation is pseudoreal leading to matrix elements 
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of the Dirac operator that are reaL Finally, for gauge group SU{Nc) in the adjoint rep- 
resentation, the gauge field is real and the matrix elements of the Dirac operator can be 
organized into real quaternions. 

In this work we want to construct a model which describes the fluctuations of the 
smallest eigenvalues as a function of the temperature. Near the theory of critical 
phenomena tells us that the fluctuations are universal with non-trivial critical exponents 
which are related to the propagation of soft modes. However, the recent work of Kocic 
and Kogut suggests that this scenario may not be valid for phase transitions involving 
bosons composed of fermions. Instead, the lowest Matsubara frequency, ttT, cures the 
infrared divergences and leads to a mean-field-like second-order phase transition. Ac- 
cording to their work on the three dimensional Gross-Neveu model, the dynamics of the 
phase transition are determined by the lowest Matsubara frequency. In this spirit, the 
only temperature effect we include in our model is that of the lowest Matsubara frequency 
between each pair of suitably chosen basis states (1 ± 75) 0„. In the sector of topological 
charge u our basis must be complemented by u unpaired basis states of the same chirality. 
Together with the symmetries mentioned above, this leads to the random matrix model, 

(2.3) 

where Vl^ is an n x [n + u) matrix. The integration over W is to be performed according 
to the Haar measure. We also include an arbitrary complex mass matrix m with mass 
eigenvalues equal to rrif. For QCD with = 2 the matrix elements of W are real 
(/5 = 1). They are complex for > 3 {l3 = 2). In each case we include fermions in the 
fundamental representation. 

The fermion determinant in ( |2.tj| ) can be written as a Grassmann integral. 

Z,iu,Nj) = I VWVrVi;eM^i:^i^'* ( ^.'f^j ^_tm^ ) V^'^ exp[-^ Tr W^W^t] 

The quenched approximation is obtained from this model as in the replica trick. We 
calculate a property for arbitrary Nf and take the limit Nf at the end of the 
calculation. In particular, quantities which are A'^y-independent are valid for Nf = as 
well. 



[2A) 
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It should be stressed that the partition function ( p.3|) represents a schematic model for 
the chiral phase transition. Although the temperature dependence of this model does not 



coincide with that of chiral perturbation theory |]TT|, it will be shown below that there is 
considerable qualitative agreement. 



3 Analytical results 



In this section we evaluate the partition function (|2.4|) using methods which are standard 
in the supersymmetric formulation of random matrix theory [12, 13] . The first step is 
to perform the average over W by performing a Gaussian integral. This leads to a four- 
fermion interaction. After averaging over the matrix elements of the Dirac operator, the 
partition function becomes 



for /? = 1, and 



Z2(z/, Nf) = J r'V^*r'^exp[--^V'{lV^mV^riV'ifc + + nifgijlliji^ 

for (3 = 2. In both cases, each of the four-fermion terms can be written as the difference of 
two squares. Each square can be linearized by the Hubbard-Stratonovitch transformation 
according to 

r 

expi-AQ"^)^ daexp{-—-iQa) . (3.3) 

J 4^1 

For j3 = 2 the partition function, expressed in terms of the two bosonic variables, can be 
combined into a single complex Nf x Nf matrix, A, resulting in 

Z2{u,Nf) = J VAVipVip* exp[ ^Ti AA^ 

+ ^{;^i,(A + m) + 4:^|,(At + m*)+^7rT(^^:^{, + ^{;^^,)] (3.4) 



For /3 = 1 six new bosonic matrix variables are required. This is related to the fact 
that, for two colors, baryons are composed of two quarks and are bosons. They can be 
combined into one antisymmetric, complex 2Nf x 2Nf matrix A resulting in the partition 
function 

xexp[^7rT(^^:V^{, + V'{:V^L)]- (3.5) 
In this case the mass matrix is an antisymmetric matrix given by 

Note also that the temperature-dependent term can be rewritten as 

Using this, the fermionic integrals can be performed, and the partition function is 
given by 

Z,{u, Nf) = J VAexp[-^TTAA^]de&\{A + m)det" ^ t/^ /'^^^ ) (3.8) 
for P = 2, and 

Z,iu,N,) = J VAeM-^TrAA^]Pf^^K-A + M)Pr ( ^'_+^* ) (3.9) 

for j3 = 1. In ( p.8| ) A is an arbitrary complex matrix whereas in ( |3.9| ) A is an arbitrary 
antisymmetric complex matrix. 

In each case the condensate is given by 

{qq) = :r\-dm\ogZ , (3.10) 

where Z is evaluated for a diagonal mass matrix with equal diagonal matrix elements. 
In the limit n ^ oo, the condensate can be determined with the aid of a saddle point 
approximation. The saddle point equations for [3 = 2 are given by 

- ——A + n{A + m)({A^ +m){A + m) +-k'^T'^) =0 . (3.11) 



An arbitrary complex matrix can be diagonalized by performing the decomposition 

A = UAV-\ (3.12) 

with all eigenvalues positive and U and V unitary matrices. We find that the solution of 
( p. lip yields U = V = 1 with eigenvalues A given by the positive roots of 

j:^X{{X + my + 7r^T^)-X-m = . (3.13) 

In the chiral limit we find a critical point at 

Tc = ^ . (3.14) 

In order to calculate the condensate, we express the derivative of the partition function 
in (3.10) in terms of an average over A, 

Below Tc the mass m can be neglected in the saddle point equation, and we find 

(gg) = S(l - vr^T^S")^/" . (3.16) 

At Tc the solution of the saddle point equation develops a non-analytic dependence on m 
resulting in the condensate 

(qq) = Stm5 . (3.17) 

Therefore, we reproduce the mean field value for the critical exponent 6 = 3. 
For (3 = 1 the saddle point equation for the A integration is 

^^\ + -{A + M)(iA^ + M){A-M) + TT^T^Y^ = . (3.18) 



2 2 

This equation can be solved by diagonalizing the complex antisymmetric matrix A as 
A = UAU, where f/ is a unitary matrix. Here, A is an antisymmetric, standard matrix 
such that Afc fc+i = — A^+i ^ = for k = 1, ■ ■ ■ , 2Nf — 1 with all other matrix elements 
equal to zero. A suitable re-definition of U can always made such that all A^ > 0. The 
condensate is calculated as in the case of (3 = 2 with the same result both below and at 
Tc. Of course, it comes as no surprise that we obtain the mean field value for the critical 
exponent in this case as well. 
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4 The phase transition 



In the remainder of this paper we describe the results of numerical investigations of the 
random matrix model (2.3). For Nj ^ 0, the determinantal weight must be included in 
the integration measure which is extremely costly. However, for Nf = the distribution 
functions are simple Gaussians, and only this case will be studied. Instead of a Gaussian 
distribution of matrix elements, we have used a rectangular distribution centered about 
zero with variance given by 

1 



where n is the size of the off-diagonal blocks. If we compare this to the variance of the 
Gaussian distribution (2.3), we can make the identification 

1 2 



6n n/3E2' ^^"^^ 

which yields = 6 for f3 = 2. According to (3.14) the critical temperature is thus 
ttTc = 1/VQ = 0.40824. General universality arguments [|l4| imply that the shape of the 



distribution of the matrix elements does not affect the properties of our random matrix 
model, and we expect that the results for the rectangular distribution will be in complete 
agreement with those for the Gaussian distribution. 

The chiral order parameter of our schematic model is the spectral density p(0). It is 
related the the chiral condensate via the Banks-Casher formula HTSl 



{qq) = ^^pio) . (4.3) 

Let us first consider the complete spectral density p(A). At T = it can be shown, using 
arguments familiar from random matrix theory |]16[, that this density has a semicircular 



shape. For non-zero temperature, it is somewhat more difficult to obtain the level density 
analytically [0. In Figs, la, lb, and Ic, we show numerical results for the complete 
spectral density for vrT = 0, vrT = 0.4, and vrT = 1.0. The deviations from a semi-circle at 
T = are not statistical fluctuations but rather well-understood finite-ra corrections ITT 



When T is large, we find that the distribution splits into two semicircles with centers at 
+ttT and — ttT. It is elementary to demonstrate that all eigenvalues come in pairs located 
symmetrically about zero. When T = 0, the distribution of the eigenvalue of smallest 
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magnitude is also known analytically. For the moment, we note only that p(A) ~ A for 
values of A less than the average value of this eigenvalue of smallest magnitude. Thus, 
some care is required in order to extract a meaningful {i.e., non-zero) value of p(0). The 
simplest method is to construct p(0) from the fraction of eigenvalues in an interval 2 AA 
centered about 0. When A A is less than the average of eigenvalue of smallest magnitude, 
this estimate of p(0) is proportional to AA. 




Figure 1: Histograms of the complete spectrum of eigenvalues obtained for n = 20 at 
temperatures of (a) nT — (top curve), (b) ttT = 0.4 (middle curve), and (c) ttT = 1.0 
(bottom curve). Each spectrum was obtained from 10^ matrices. 

For somewhat larger values of AA there is a relatively rapid saturation of p(0) to its desired 
asymptotic value. In principle,it is better to obtain p(0) by performing a smearing of the 
full spectrum with a Gaussian whose width is roughly comparable to the spacing between 
adjacent eigenvalues. This avoids the small- amphtude oscillations associated with the 
first approach. In practice, there is little difference between the two methods. 

According to the mean field argument presented in section 3, our model shows a 
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second-order phase transition with 

p(0, T) ~ - r2 . (4.4) 

In order to account for finite n effects approximately, we convolute this expression with a 
Gaussian 

p,(0, T) ~ dxJr^ - x2 exp[-(x - Tf/a^] . (4.5) 

While we offer no analytic justification for this form, we note that it provides an excellent 
fit to the results of simulations. A best fit of this expression to our results for = 20 is 
shown in Fig. 2. The spectral density was obtained from a bin size of 0.02. In this case, 
the fitted values for Tc and a are 0.389 and 0.074, respectively. Using data for n > 50 
we performed an extrapolation to n ^ oo. We find a critical temperature of Tc = 0.3967 
which is in good agrees with the theoretical value of vrTc = 1/ for our present numerical 
parameters. 




0.00 0,20 0.40 0,60 

T 

Figure 2: The (qq) condensate as calculated from (4.3) for n = 20. Each point was 
obtained on the basis of 2 x lO'' matrices, and p(0) was obtained by counting eigenvalues 
in the range < A < 0.02. The solid curve is a fit according to (4.5) using the parameters 
Tc = 0.3890 and a = 0.7414. 

Recent lattice QCD calculations have studied the condensate as a function of the 
so-called valence mass. In these lattice simulations, the sea-quark mass was much larger 
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than the smallest eigenvalue so that the problem is effectively equivalent to taking Nf = 0. 
It is possible to simulate these calculations within the framework of the present model by 
defining 



1 " 

E 



m 



(4.6) 



Fig. 3 shows a log-log plot of the resulting values for (ijjip) as a function of m for n = 20 
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Figure 3: The condensate as a function of the so-called valence quark mass calculated 
according to (4.6) at ttT = 0, ttT = 0.4 and ttT = 1.0. The calculations were performed 
for n = 20 and 2 x 10^ matrices were used for each temperature. 

Our results are qualitatively similar to the lattice calculations. The m —>■ limit is clearly 
ambiguous. These results show that the spectral density yields a much more accurate 
determination of the chiral condensate Ijl^ . 



5 Dynamics of the smallest eigenvalue 

In this section we study the distribution of the eigenvalue of smallest magnitude, denoted 
by Ai, as a function of temperature. This distribution is known analytically for a few 
special cases of the Laguerre ensemble. The present case, defined by the random matrix 
model (2.3) for Nf = 0, at T = happens to be one of them. With our choice of 
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parameters, the distribution is given by [|0 



/(A0~A.exp('-(^«^2)') . (5.1) 

This T = result is exact even for finite n. (The factor of Ai is readily understood. 
For every eigenvalue +A there is always a corresponding eigenvalue of —A. This factor is 
simply a consequence of level repulsion.) In order to be able to describe the temperature 
dependence, we introduce a more general distribution 

/(Ai) = A, exp (-^^^^) . (5.2) 

At zero temperature xq = and a = 2/n /5 S. The average value of the smallest eigenvalue 
is given as Ai = ^ jnfi S. While this expression is not, in general, exact, it is serviceable. 




Figure 4: The distribution of the eigenvalue of smallest magnitude determined from 
2 X 10^ trials with n = 20 for a bin size of 0.02 (points). The histogram was obtained 
from (5.2) with a;o = 0. The value of o" = 0.020685 was determined from the average value 
of Ai obtained from the simulations. This fitted value is in excellent agreement with the 
analytic value of 1jn(5Yj = 0.02041 . . .. The form of (5.2) provides an excellent fit to the 
simulation data with = 0.8. 

In Fig. 4 we show a histogram of the distribution of the smallest eigenvalue at T = 
for the case n = 20. The results are in perfect agreement with (4.1). The temperature 
dependence of the parameters of the distribution is shown in Fig. 5, where we show the 



12 



ratio Ai/AAi as a function of the temperature for n — 10. When xo — this ratio is 

equal to [7r/(4 — tt)]^/^ = 1.913 Surprisingly, we find that that this ratio is constant 

for T < Tc. In terms of xq, this implies that xq is strictly zero for T < Tc. Inspection 
of the corresponding distribution indicates that (5.1) remains quantitatively valid in this 
region. Above this temperature (and ignoring some finite n effects near threshold), Xq 
grows linearly with T. This behavior is expected since the entire distribution moves 
hnearly with T for sufficiently large T. What is surprising is that xq goes to zero below 
some finite T and that the effect is so pronounced for a matrix of such small dimension. 




Figure 5: The ratio of the smallest eigenvalue to its rms deviation, Ai/AAi, as a function 
of temperature for n = 10. Each point represents 2 x 10^ matrices. The average value 
of this ratio for T < Tc is 1.9107 which is in good agreement with the expected value of 
(7r/(4-7r))V2 = 1.913.... 



The behavior of Xo{T) versus T (shown in Fig. 6) allows us to extract from x{Tc) = 0. 
In practice we use a linear extrapolation of our data for T > Tc. Our results for different 
size matrices can be summarized by the expression 



ri + 2.51316 
n + 4.10346 ' 



(5.3) 



which yields an asymptotic result of ttTc = 0.40859, which is very close to the theoretical 
result of l/-\/6 = 0.40824. As indicated, the form of (5.2) is not exact. Similar results can 
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be obtained in a model-independent fashion by straight-line extrapolation of the observed 
values of Ai/AAi as a function of the temperature. 
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Figure 6: The parameter xq of (5.2) as a function of T for n = 10. Each point represents 
2 X 10"^ matrices. The parameters Xq and a were determined from the ratio Ai/AAi (shown 
in Fig. 5) and the corresponding values of Ai. Linear extrapolation suggests ttTc = 0.3625 
for this case. 

Finally, wc study the scaling of the smallest eigenvalue with N. For T < Tc, the Banks- 
Casher formula suggests that Amin ~ 1/n. Indeed, this is what is found numerically both 
for the expectation value and the variance of Amin- 

At T — Tc we have that E(m) ~ m^^^. This leads to the eigenvalue density (for A — > 0) 



Smaller masses suppress the eigenvalue density near zero, so we must have cti > 0. For 
a fixed mass, the eigenvalue density should not diverge, so we must also have a2 > 0. It 
is possible that ai and 02 depend on Nf, but we were not able to investigate this point 
within the present framework. The scaling behavior of the smallest eigenvalue is obtained 
from Jq'^'" p{X)dX ^ l/N which leads to 



p(A)~m"iA"2 with cti + ^2 



1 

~5' 



(5.4) 




(5.5) 
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For Nf = 0, the eigenvalue density cannot depend on m and we thus have 

a2 = - and ai = for Nf = . (5.6) 

This yields (Amin) ~ A^"^''^, which has been verified numerically to very high accuracy. It 
was also shown that the variance of Amin has the same scaling behavior. 

For T > Tc the average position of the smallest eigenvalue departs from zero and scales 
like N°. We expect ^ that its variance well have the scaling behavior which is 

typical for eigenvalues near the edge of a semicircle. Again, this is in perfect agreement 
with our numerical simulations. 



6 Conclusions 

In this paper we have studied a random matrix model possessing the global symmetries 
of the QCD-action and the temperature dependence suggested by the form of the lowest 
Matsubara frequency. At T = this model is completely soluble; it reduces to what 
is known in the mathematical literature as the Laguerre ensemble. For nonzero tem- 
peratures, we have succeeded in extracting some interesting properties analytically. In 
particular, we have shown that the model undergoes a second-order phase transition. Two 
critical exponents were obtained. 

6 = 3. (6.1) 

The lattice result [jl| for QCD with two light fiavors for 1/(36 is 0.77±0.14, which is closer 
to these mean field values than to the results for either the 0(4) or 0(2) Heisenberg spin 
models. For three or more light fiavors, QCD shows a first-order chiral phase transition, 
whereas for one fiavor there is no transition at all. Our model does not contain such 
fiavor dependence: It has a second-order phase transition for any number of fiavors. If 
our model has anything to say about QCD, it is for QCD with two light fiavors. 

Numerically, we found one surprising result: The distribution of the smallest eigenvalue 
below Tc is (numerically) equivalent to the distribution obtained for T = 0. For T > T^, 
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the centroid of the Gaussian distribution grows hnearly with T. The behavior at T = 
agrees well with known analytical results. 

Our results suggest an alternative method for obtaining the critical temperature; 
namely from the distribution of the eigenvalue of smallest magnitude. It would be in- 
teresting to study the dynamics of the smallest eigenvalues in lattice QCD as well. 



Acknowledgements 

The reported work was partially supported by the US DOE grant DE-FG-88ER40388. 
We would hke to thank H.A. Weidenmiiller and T. Wettig for useful discussions and for 
communicating their results on a related model prior to publication. 



16 



References 

[I] F. Karsch and E. Laermann, Susceptibilities, the specific heat and a cumulant in two- 
flavor QCD, Bielefeld preprint BI-TP 94/29. 

[2] S. Chandrasekharan, in Continuous Advances in QCD, eds. A. Smilga, World Scien- 
tific, Singapore 1994. 

[3] R.D. Pisarski and F. Wilczek, Pliys. Rev. D29 (1984) 338. 

[4] K. Rajagopal and F. Wilczek, Nucl. Pliys. B399 (1993) 395. 

[5] A. Kocic and J. Kogut, Pliys. Rev. Lett. 74 (1995) 3109; Phase transitions at finite 
temperature and dimensional reduction for fermions and boson. University of Illinois 
preprint, 1995 ( |hep-lat/9507CTT^ ). 

[6] S.J. Hands and M. Teper, Nucl. Pliys. B347 (1990) 819. 

[7] T. Kalkreuter, Pliys. Rev. D, (in press), |hep-lat/94080T^ ; Pliys. Lett. B276 (1992) 
485; Phys. Rev. D48 (1993) 1926. 

[8] J. Verbaarschot, Phys. Rev. Lett. 72 (1994) 2531; Pliys. Lett. B329 (1994) 351; Nucl. 
Phys. B426 [FS] (1994) 559; Nucl. Phys. B427 (1994) 434. 

[9] A. Halasz and J.J.M. Verbaarschot, Effective Lagrangians and chiral random matrix 
theory, Phys. Rev. D (1995) (in press). 

[10] H. Leutwyler and A. Smilga, Phys. Rev. D46 (1992) 5607. 

[II] J. Gasser and H. Leutwyler, Phys. Lett. B188 (1987) 477. 
[12] K.B. Efetov, Adv. Mod. Phys. 32 (1983)53. 

[13] J. Verbaarschot, H. Weidenmiiller, and M. Zirnbauer, Phys. Rep. 129 (1985) 367. 

[14] E. Brezin and A. Zee, Nucl. Phys. B402 (1993) 613; C. Beenhakker, Nucl. Phys. 
B422 (1994) 515; G. Hackenbroich and H. Weidenmiiller, Phys. Rev. Lett. (1995). 

[15] T. Banks and A. Casher, Nucl. Phys. B169 (1980) 103. 

17 



[16] M. Mehta, Random Matrices, Academic Press, San Diego, 1991 

[17] J.J.M. Verbaarschot and I. Zahed, Phys. Rev. Lett. 70 (1993) 3852. 

[18] A.D. Jackson, M. Sener and J.J.M. Verbaarschot, in preparation. 

[19] J.J.M. Verbaarschot, Universal scaling of the valence quark mass dependence of the 
chiral condensate, preprint SUNY-NTG-95/32. 

[20] P.J. Forrester, Nucl. Phys. B 402[FS] (1993) 709. 

[21] J.J.M. Verbaarschot and M.R. Zirnbauer, Ann. Phys. 158 (1984) 78. 



18 



